home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / network / blast2 / client / blstout2.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-07-05  |  42.2 KB  |  1,513 lines  |  [TEXT/R*ch]
  1. /*
  2. * ===========================================================================
  3. *
  4. *                            PUBLIC DOMAIN NOTICE                          
  5. *               National Center for Biotechnology Information
  6. *                                                                          
  7. *  This software/database is a "United States Government Work" under the   
  8. *  terms of the United States Copyright Act.  It was written as part of    
  9. *  the author's official duties as a United States Government employee and 
  10. *  thus cannot be copyrighted.  This software/database is freely available 
  11. *  to the public for use. The National Library of Medicine and the U.S.    
  12. *  Government have not placed any restriction on its use or reproduction.  
  13. *                                                                          
  14. *  Although all reasonable efforts have been taken to ensure the accuracy  
  15. *  and reliability of the software and data, the NLM and the U.S.          
  16. *  Government do not and cannot warrant the performance or results that    
  17. *  may be obtained by using this software or data. The NLM and the U.S.    
  18. *  Government disclaim all warranties, express or implied, including       
  19. *  warranties of performance, merchantability or fitness for any particular
  20. *  purpose.                                                                
  21. *                                                                          
  22. *  Please cite the author in any work or product based on this material.   
  23. *
  24. * ===========================================================================
  25. *
  26. * File Name: blstout2.c
  27. *
  28. * Author:  Roman L. Tatusov, Warren Gish, Jonathan Epstein, Tom Madden, Yuri Sadykov
  29. *
  30. * Version Creation Date:   06/16/95
  31. *
  32. * $Revision: 4.0 $
  33. *
  34. * File Description: 
  35. *       Creating BLAST report
  36. *       Warren Gish's utilities needed for producing BLAST output
  37. *
  38. * Modifications:  
  39. * --------------------------------------------------------------------------
  40. * Date     Name        Description of modification
  41. * -------  ----------  -----------------------------------------------------
  42. *
  43. * ==========================================================================
  44. *
  45. *
  46. * RCS Modification History:
  47. * $Log: blstout2.c,v $
  48.  * Revision 4.0  1995/07/26  13:55:34  ostell
  49.  * force revision to 4.0
  50.  *
  51.  * Revision 1.4  1995/07/19  15:00:44  madden
  52.  * Changed some spacings in TraditionalHistBlastOutput.
  53.  *
  54.  * Revision 1.3  1995/07/18  21:06:51  madden
  55.  * changed p-value etc. to p_value etc.
  56.  *
  57.  * Revision 1.2  1995/06/22  17:03:44  madden
  58.  * HitDataPtr replaced by BLAST0ResultPtr.
  59.  *
  60.  * Revision 1.1  1995/06/16  11:26:47  epstein
  61.  * Initial revision
  62.  *
  63. */
  64.  
  65. #include <ncbi.h>
  66. #include <objseq.h>
  67. #include <objsset.h>
  68. #include <sequtil.h>
  69. #include <seqport.h>
  70. #define NLM_GENERATED_CODE_PROTO
  71. #include <netblap2.h>
  72. #include <objblst2.h>
  73.  
  74. #define BLASTCLI_MATRIX_SIZE 26
  75.  
  76. static char ralpha[256];
  77. static Int4 matr[BLASTCLI_MATRIX_SIZE+1][BLASTCLI_MATRIX_SIZE+1];
  78. static int prots;
  79. double la2;
  80.  
  81. /*************************************************************
  82. *                                                             *
  83. *  Warren Gish's utilities needed for producing BLAST output  *
  84. *                                                             *
  85.  *************************************************************/
  86.  
  87. /* wrap -- wordwrap lines of output */
  88.  
  89. static void    doindent PROTO((FILE *fp, int ncols));
  90. static CharPtr    gish_str_strand PROTO ((Uint2 strand));
  91. static Int2 format_an_id PROTO ((CharPtr id, ValNodePtr vnp, Int2 max_id_length));
  92. static Int2 get_max_ID_length PROTO ((BLAST0HitListPtr hlp));
  93. static void get_frame PROTO ((BLAST0SeqIntervalPtr loc, CharPtr array, Int4 total_length));
  94.  
  95. static int _cdecl
  96. wrap(fp, title, s, slen, linelen, indent)
  97.     FILE    *fp;    /* the output stream */
  98.     char    *title;    /* title string */
  99.     char    *s;    /* pointer to null-terminated string for output */
  100.     int        slen;    /* strlen(s), or -1 */
  101.     int        linelen,    /* max. length of an output line before each '\n' */
  102.             indent;        /* no. of columns to indent any continuation lines */
  103. {
  104.     register char    *savep, *savep2;
  105.     register char    *cp, ch;
  106.     int        outlen, len, olinelen;
  107.     int        titlelen;
  108.     Boolean    once = TRUE;
  109.  
  110.     if (title != NULL) {
  111.         cp = title;
  112.         len = 0;
  113.         while ((ch = *cp++) != NULLB) {
  114.             fputc(ch, fp);
  115.             if (ch == '\n')
  116.                 len = cp - title;
  117.         }
  118.         titlelen = (cp - title) - len - 1;
  119.     }
  120.     else
  121.         titlelen = 0;
  122.  
  123.     if (slen < 0)
  124.         slen = strlen(s);
  125.  
  126.     if (indent >= linelen) {
  127.         indent = MIN(1 + titlelen, linelen-5);
  128.         indent = MAX(indent, 0);
  129.     }
  130.     olinelen = linelen;
  131.     linelen -= titlelen;
  132.     for (cp = s; cp < s + slen;) {
  133.         /* Skip leading white space */
  134.         for (;; ++cp) {
  135.             if ((ch = *cp) == NULLB)
  136.                 goto Done;
  137.             if (!IS_WHITESP(ch))
  138.                 break;
  139.         }
  140.  
  141.         outlen = cp - s;
  142.         if (slen - outlen <= linelen) {
  143.             /* Remainder is short enough to fit on one line */
  144.             if (!once)
  145.                 doindent(fp, indent);
  146.             fprintf(fp, "%s\n", cp);
  147.             break;
  148.         }
  149.         else {
  150.             savep2 = cp + linelen;
  151.             ch = *savep2;
  152.             if (IS_WHITESP(ch)) {
  153. Phase2:
  154.                 for (savep = savep2; savep >= cp; --savep)
  155.                     if (!IS_WHITESP(*savep)) {
  156.                         ++savep;
  157.                         break;
  158.                     }
  159.             }
  160.             else {
  161.                 for (savep = savep2; savep >= cp; --savep)
  162.                     if (IS_WHITESP(*savep)) {
  163.                         savep2 = savep;
  164.                         goto Phase2;
  165.                     }
  166.                 /* a _very_ long word here! */
  167.                 savep = savep2;
  168.             }
  169.         }
  170.         if (!once)
  171.             doindent(fp, indent);
  172.         once = FALSE;
  173.         while (cp < savep)
  174.             fputc(*cp++, fp);
  175.         putc('\n', fp);
  176.         cp = savep2;
  177.         linelen = olinelen - indent;
  178.     }
  179.  
  180. Done:
  181.     return ferror(fp);
  182. }
  183.  
  184. static void
  185. doindent(fp, ncols)
  186.     register FILE    *fp;
  187.     register int    ncols;
  188. {
  189.     while (ncols-- > 0)
  190.         fputc(' ', fp);
  191. }
  192.  
  193. /**** MISC  ****/
  194.  
  195. static CharPtr
  196. gish_str_strand(Uint2 strand)
  197. {
  198.     switch (strand)
  199.     {
  200.     case BLAST0_Seq_interval_strand_plus:
  201.     case BLAST0_Seq_interval_strand_plus_rf:
  202.         return "Plus"; /* plus */
  203.     case BLAST0_Seq_interval_strand_minus:
  204.     case BLAST0_Seq_interval_strand_minus_rf:
  205.         return "Minus"; /* minus */
  206.     case BLAST0_Seq_interval_strand_both:
  207.     default: /* unknown, or unable to convert */
  208.         return "Undefined";
  209.     }
  210. }
  211.  
  212. static BLAST0ResponsePtr _find(BLAST0ResponsePtr p, Uint1 v)
  213. {
  214.     BLAST0ResponsePtr b;
  215.  
  216.     for (b = p; b != NULL && b->choice != v; b = b->next) 
  217.         ;
  218.     return b;
  219. }
  220.  
  221. static Pointer find(BLAST0ResponsePtr p, Uint1 v)
  222. {
  223.     BLAST0ResponsePtr b;
  224.  
  225.     for (b = p; b != NULL && b->choice != v; b = b->next) 
  226.         ;
  227.  
  228.     return (b == NULL) ? NULL : b->data.ptrvalue;
  229. }
  230.  
  231. static CharPtr
  232. print_double(CharPtr pCh, double dX, int iWidth, int iPrecision)
  233. {
  234.    double dY;
  235.    int    i;
  236.    
  237.    dY = log(dX * 1.000000001) / NCBIMATH_LN10;
  238.    i = dY;
  239.    dY = Nlm_Powi(10., i - iPrecision + 1);
  240.    dX = Nlm_Nint(dX / dY) * dY;
  241.    dY = Nlm_Powi(10., iPrecision - 1) - 1.e-7;
  242.    if (dX >= dY)
  243.       sprintf(pCh, "%*.0lf", iWidth, dX);
  244.    else
  245.       sprintf(pCh, "%*.*lf", iWidth, iPrecision - i - 1, dX);
  246.    
  247.    return pCh;
  248. }
  249.  
  250. /*
  251. The following function serves a very important purpose.  It translates
  252. the different alphabets into ncbistdaa.  This is necessary for 
  253. calculations involving the BLOSUM/PAM matrics.
  254. */
  255.  
  256. static void get_ralpha(char *ralpha, Uint1 code)
  257. {
  258.     SeqCodeTablePtr sctp;
  259.     int index_i, index_j;
  260.  
  261.     sctp = SeqCodeTableFind(code);
  262.     for (index_i = 0; index_i < 256; index_i++) {
  263.         index_j = GetResidueForSymbol(sctp, index_i);
  264.         if (index_j == 255) 
  265.             index_j = 0;
  266.         ralpha[index_i] = index_j;
  267.     }
  268. }
  269.  
  270. /***************************************************************************
  271. *
  272. *    Format the top part of the "traditional" BLAST output.
  273. *    This includes the name of the program, the version, the compile
  274. *    date, the reference down through "Searching...":
  275. *
  276. *BLASTN 1.4.7MP [16-Oct-94] [Build 16:01:48 Dec 16 1994]
  277. *
  278. *Reference:  Altschul, Stephen F., Warren Gish, Webb Miller, Eugene W. Myers,
  279. *and David J. Lipman (1990).  Basic local alignment search tool.  J. Mol. Biol.
  280. *215:403-10.
  281. *
  282. *Notice:  this program and its default parameter settings are optimized to find
  283. *nearly identical sequences rapidly.  To identify weak similarities encoded in
  284. *nucleic acid, use BLASTX, TBLASTN or TBLASTX.
  285. *
  286. *Query=  195747|Mitochondrion Rupicapra rupi
  287. *        (646 letters)
  288. *
  289. *Database:  Non-redundant PDB+GBupdate+GenBank+EMBLupdate+EMBL
  290. *           248,447 sequences; 236,652,708 total letters.
  291. *Searching..................................................done
  292. *
  293. *    The matrix used for the BLAST run is also read in.
  294. *    
  295. **************************************************************************/
  296. static Boolean
  297. TraditionalHeadBlastOutput(BLAST0ResponsePtr blresp, Int4Ptr query_length, CharPtr program, FILE *fp)
  298. {
  299.     BLAST0ResponsePtr a;
  300.     Scores_scaled_intsPtr sco;
  301.        ValNodePtr   ints, sc;
  302.     BLAST0MatrixPtr ma;
  303.     BLAST0KABlkPtr ka;
  304.     BLAST0StatusPtr status;
  305.     BLAST0PrefacePtr pre;
  306.     BLAST0SequencePtr qu;
  307.     BLAST0DbDescPtr db;
  308.     BLAST0JobDescPtr jobd;
  309.     Int2 index1, index2, length;
  310.     Int4 progr;
  311.     Char buf[30];
  312.     CharPtr string;
  313.     ValNodePtr node;
  314.     CharPtr defl;
  315.     
  316.     if ((pre = find(blresp, BLAST0Response_preface)) != NULL) {
  317.         prots = (pre->program[5] != 'N');
  318.         fprintf(fp, "%s %s [%s]", pre->program, pre->version, pre->dev_date);
  319.         if (pre->bld_date != NULL) {
  320.             fprintf(fp, " [Build %s]", pre->bld_date);
  321.         }
  322.         fprintf(fp, "\n");
  323.         for (node = pre->cit; node != NULL; node = node->next) {
  324.             string = (CharPtr) node->data.ptrvalue;
  325.             if (node == pre->cit)
  326.                 wrap(fp, "\nReference:  ", string, strlen(string), 79, 0);
  327.             else
  328.                 wrap(fp, "", string, strlen(string), 79, 0);
  329.         }
  330.         for (node = pre->notice; node != NULL; node = node->next) {
  331.             string = (CharPtr) node->data.ptrvalue;
  332.             wrap(fp, "\nNotice:  ", string, strlen(string), 79, 0);
  333.         }
  334.     }
  335.     if ((qu = find(blresp, BLAST0Response_query)) != NULL) {
  336.         fprintf(fp, "\nQuery=  ");
  337.         if (qu != NULL && qu->desc != NULL &&
  338.             qu->desc->defline != NULL) {
  339.             node = qu->desc->id;
  340.             while (node)
  341.             {
  342.             if (node->choice == BLAST0SeqId_giid)
  343.             {
  344.                 sprintf(buf, "gi|%ld|", (long) node->data.intvalue);
  345.                 string = &buf[0];
  346.             }
  347.             else if (node->choice == BLAST0SeqId_textid)
  348.             {
  349.                     string = node->data.ptrvalue;
  350.                     if (strncmp(string, "lcl|", 4) == 0)
  351.                 {
  352.                     sprintf(buf, "%s ", string+4);
  353.                     string = &buf[0];
  354.                 }
  355.             }
  356.                 fprintf(fp, "%s", string);
  357.             node = node->next;
  358.             }
  359.             defl = qu->desc->defline;
  360.             wrap(fp, " ", defl, strlen(defl), 79, 8);
  361.         } else {
  362.             fprintf(fp, "Unknown\n");
  363.         }
  364.         *query_length = qu->length;
  365.         fprintf(fp, "        (%s letters)\n", Nlm_Ultostr(qu->length, 1));
  366.     }
  367.  
  368.     if (StringCmp(program, "tblastx") == 0)
  369.     {
  370.         fprintf(fp, "\n  Translating both strands of query sequence in all 6 reading frames\n");
  371.     }
  372.  
  373.     if ((db = find(blresp, BLAST0Response_dbdesc)) != NULL) {
  374.         string = db->def;
  375.         wrap(fp, "\nDatabase:  ", string, strlen(string), 79, 11);
  376.         fprintf(fp, "           %s sequences; %s total letters.\n",
  377.             Nlm_Ultostr((long)(db->count), 1), Nlm_Ultostr(db->totlen, 1));
  378.     }
  379.  
  380.     a = blresp;
  381.     while ((a = _find(a, BLAST0Response_job_start)) != NULL) {
  382.         jobd = a->data.ptrvalue;
  383.         if (jobd->jid != BLAST0_Job_desc_jid_search) {
  384.             a = a->next;
  385.             continue;
  386.         }
  387.         fprintf(fp, "%s", jobd->desc);
  388.         fflush(fp);
  389.         progr = BLAST0Response_job_progress;
  390.         for (a = a->next; a != NULL && a->choice == progr; a = a->next) {
  391.             fprintf(fp, ".");
  392.             fflush(fp);
  393.         }
  394.         if (a != NULL && a->choice == BLAST0Response_job_done) {
  395.             fprintf(fp, "done\n");
  396.         } else {
  397.             fprintf(fp, "Interrupted!!!\n");
  398.         }
  399.         fflush(fp);
  400.     }
  401.     
  402.     if ((status = find(blresp, BLAST0Response_status)) != NULL) {
  403.         if (status->code != 0) {
  404.             fprintf(fp, "ERROR: %ld\n", (long) status->code);
  405.             fprintf(fp, "reason: %s\n", status->reason);
  406.         }
  407.     }
  408.     fflush(fp);
  409.     if ((ma = find(blresp, BLAST0Response_matrix)) != NULL) {
  410.         sc = ma->Scores_scores;
  411.         sco = sc->data.ptrvalue;
  412.         ints = sco->ints;
  413.         if (ints)
  414.         {
  415.            for (index1=0; index1<BLASTCLI_MATRIX_SIZE; index1++)
  416.                for (index2=0; index2<BLASTCLI_MATRIX_SIZE; index2++)
  417.                {
  418.                   matr[index1][index2] = ints->data.intvalue;
  419.                if (ints->next == NULL)
  420.                {    /* Get out of both loops! */
  421.                 index1 = BLASTCLI_MATRIX_SIZE;
  422.                 index2 = BLASTCLI_MATRIX_SIZE;
  423.                 break;
  424.                }
  425.                ints = ints->next;
  426.                }
  427.         }
  428.         get_ralpha(ralpha, Seq_code_ncbistdaa);
  429.     }
  430.     if ((ka = find(blresp, BLAST0Response_kablk)) != NULL) {
  431.         la2 = ka->lambda / log(2.);  /* lambda / ln2  */
  432.     }
  433.     return TRUE;
  434. }
  435.  
  436. #define PAGE_W 80
  437. #define EXPECT_W 6
  438. #define EXPECT_PRCSN 3
  439. #define SEPARATOR_SYMBOL '|'
  440. #define DRAWING_SYMBOL '='
  441. #define SMALLDRAWING_SYMBOL ':'
  442. /**************************************************************************
  443.  *    The function formats histgram data and is part of the "traditional"
  444.  * BLAST output.
  445.  * 
  446.  * Note:
  447.  *    If pointer on the histogram structure == NULL it doesn't seem as
  448.  * error, simply it means no histogram data have been requested. In this
  449.  * case corresponding message is printed to output.
  450.  * 
  451.  **************************************************************************/
  452. static Boolean
  453. TraditionalHistBlastOutput(BLAST0ResultPtr hdp, FILE *fp)
  454. {
  455.    BLAST0HistogramPtr pHist = hdp->hist;
  456.    Int4 nBars;     /* Number of histogram bars */
  457.    Int4 i;
  458.    Int4 nMaxObserved;
  459.    Int4 nMaxHist;
  460.    Int4 nHistValue;
  461.    Int2 nObsWidth; /* Width (in symbols) of field wich represents observed value */
  462.    Int2 nHistWidth;/* Width (in symbols) of field wich represents histogram value */
  463.    Int2 nCols;     /* Width of space (in symbols) where histogram is painted */
  464.    Int2 nPerCol;   /* Number of sequences per unit */
  465.    Int2 nUnits;    /* Number units needed to represent histogram bar */
  466.    Int2 j;
  467.    FloatLo fDelta;
  468.    Char aCh[20];
  469.    Boolean bNeedCheck = TRUE;
  470.    BLAST0HistogramBarPtr pBar;
  471.       
  472.    
  473.    if (pHist == NULL)
  474.       return TRUE;
  475.    
  476.    
  477.    putc('\n', fp);
  478.    putc('\n', fp);
  479.    fprintf(fp, "     Observed Numbers of Database Sequences Satisfying\n");
  480.    fprintf(fp, "    Various EXPECTation Thresholds (E parameter values)\n");
  481.    putc('\n', fp);
  482.  
  483.    nBars  = pHist->nbars;
  484.    
  485.    if (nBars == 0) {
  486.       fprintf(fp, "\n\n***  histogram slots are all empty ***\n\n");
  487.       return TRUE;
  488.    }
  489.       
  490.    /* Find out max observed and histogram data values
  491.     */
  492.    for (i = 0, nMaxObserved = 0, nMaxHist = 0, pBar = pHist->bar;
  493.     i < nBars; i++) {
  494.       nMaxObserved = MAX(nMaxObserved, pBar->n);
  495.        
  496.       if (i == nBars - 1) /* the last one */
  497.      nMaxHist = MAX(nMaxHist, pBar->n - pHist->base);
  498.       else
  499.      nMaxHist = MAX(nMaxHist, pBar->n - pBar->next->n);
  500.        
  501.       pBar = pBar->next;
  502.    }
  503.    
  504.    nObsWidth = Nlm_Ulwidth(nMaxObserved, 0);
  505.    nObsWidth = MAX(nObsWidth, 2);
  506.    nHistWidth = Nlm_Ulwidth(nMaxHist, 0);
  507.    nHistWidth = MAX(nHistWidth, 2);
  508.    
  509.    nCols = PAGE_W - (EXPECT_W+1 + nObsWidth+1 + nHistWidth+1 + 2);
  510.    fDelta = (float)nMaxHist / (float)nCols;
  511.    nPerCol = ceil(fDelta);
  512.    nPerCol = MAX(nPerCol, 1);
  513.    
  514.    fprintf(fp, "        Histogram units:      %c %d Sequence%s",
  515.        DRAWING_SYMBOL, nPerCol, (nPerCol == 1 ? "" : "s"));
  516.  
  517.     if (fDelta > 1.0)
  518.       fprintf(fp, "     %c less than %d sequences\n", SMALLDRAWING_SYMBOL, nPerCol);
  519.     else
  520.       putc('\n', fp);
  521.     
  522.    putc('\n', fp);
  523.    fprintf(fp, " EXPECTation Threshold\n");
  524.    fprintf(fp, " (E parameter)\n");
  525.    fprintf(fp, "    |\n");
  526.    fprintf(fp, "    V   Observed Counts-->\n");
  527.    
  528.    /* Draw the histogram
  529.     */
  530.    for (i = 0, pBar = pHist->bar; i < nBars; i++) 
  531.    {
  532.       if ( bNeedCheck && (pBar->x <= pHist->expect + 1e-6)) {
  533.      /* This line must appear only once
  534.       */
  535.      fprintf(fp,
  536.          " >>>>>>>>>>>>>>>>>>>>>  Expect = %#0.3lg, Observed = %lu  <<<<<<<<<<<<<<<<<\n",
  537.          pHist->expect, pHist->observed);
  538.      bNeedCheck = FALSE;
  539.       }
  540.  
  541.       nHistValue = (i < nBars - 1) ? (pBar->n - pBar->next->n) : (pBar->n - pHist->base);
  542.       fprintf(fp, " %s %*lu %*lu %c", print_double(aCh, pBar->x, 6, 3),
  543.           nObsWidth, pBar->n, nHistWidth, nHistValue, SEPARATOR_SYMBOL);
  544.  
  545.       /* Draw long line
  546.        */
  547.       if (fDelta <= 1.00) {
  548.       /* Don't stretch on whole width
  549.        */
  550.       for (j = 0; j < nHistValue; j++)
  551.         putc(DRAWING_SYMBOL, fp);
  552.       } else {
  553.       /* Stretch on whole width
  554.        */
  555.       nUnits = (FloatLo) nHistValue / fDelta;
  556.       for (j = 0; j < nUnits; j++)
  557.         putc(DRAWING_SYMBOL, fp);
  558.       }
  559.           
  560.       /* Draw tiny line if needed
  561.        */
  562.       if (j == 0 && nHistValue > 0 && fDelta > 1.00)
  563.     putc(SMALLDRAWING_SYMBOL, fp);
  564.       
  565.       putc('\n', fp);
  566.       
  567.       pBar = pBar->next;
  568.    }
  569.    
  570.    return TRUE;
  571. }
  572.  
  573. #define DEF_CUTOFF 43
  574. #define DEF_AND_ID_CUTOFF 59
  575. #define BLASTCLI_ID_MAX 25
  576.  
  577. /**************************************************************************
  578. *    
  579. *    This function formats the One-line summaries of the traditional
  580. *    BLAST output.  This includes the High-Scoring Segment Pair (HSP) 
  581. *    id, the defline of the HSP, the high score of that HSP, the P 
  582. *    value, and the number of HSP's found.  If a translation 
  583. *    was performed on the database, the frame is also printed.  An
  584. *    example is:
  585. *
  586. *                                                                     Smallest
  587. *                                                                       Sum
  588. *                                                              High  Probability
  589. *Sequences producing High-scoring Segment Pairs:              Score  P(N)      N
  590. *
  591. *dbj|D32197|RURMTCB27 Mitochondrion Rupicapra rupicapra ge...  3230  7.0e-264  1
  592. *dbj|D32191|CCRMTCB21 Mitochondrion Capricornis crispus ge...  2609  1.8e-212  1
  593. *gb|U17862|OMU17862   Ovibos moschatus moschatus cytochrom...  2600  1.0e-211  1
  594. *dbj|D32195|CCRMTCB25 Mitochondrion Capricornis sumatrensi...  2582  5.6e-210  1
  595. *gb|U17861|NCU17861   Nemorhaedus caudatus cytochrome b ge...  2573  1.8e-209  1
  596. *
  597. *    and so on.
  598. ******************************************************************************/
  599.  
  600. static Boolean
  601. TraditionalTopBlastOutput(BLAST0ResultPtr hdp, Int4 query_length, CharPtr program, FILE *fp)
  602. {
  603.     BLAST0HitListPtr hlp;
  604.     BLAST0SegmentPtr segp;
  605.     /* Int4 dim; */
  606.     BLAST0SeqIntervalPtr hit;
  607.     BLAST0SeqDescPtr sdp, sdp1;
  608.     CharPtr defline=NULL, ptr;
  609.     Char frame_array[3], buffer[DEF_CUTOFF+5];
  610.     ScorePtr score; 
  611.     CharPtr scoreDescr;
  612.     BLAST0HSPPtr hsp;
  613.     Int4 highScore, string_length;
  614.     FloatHi poissonScore;
  615.     Int2 index, max_id_length, max_def_length, length, length_left, old_length, total_length;
  616.     Int4 nScore;
  617.     Int4 TemphighScore;
  618.     FloatHi TemppoissonScore;
  619.     Int4 TempnScore;
  620.     CharPtr str1;
  621.     Char str2[DEF_CUTOFF+4];
  622.     Char id[BLASTCLI_ID_MAX+1];
  623.  
  624.     fprintf(fp, "\n%-69sSmallest\n%-71sSum\n", "", "");
  625.     if (StringCmp("blastn", program) == 0 || 
  626.     StringCmp("blastp", program) == 0)
  627.     {
  628.         fprintf(fp, "%-62sHigh  Probability\n", "");
  629.         fprintf(fp, "Sequences producing High-scoring Segment Pairs:%-14sScore  P(N)      N\n\n", "");
  630.     }
  631.     else
  632.     {
  633.         fprintf(fp, "%-53sReading  High  Probability\n", "");
  634.         fprintf(fp, "Sequences producing High-scoring Segment Pairs:%-8sFrame Score  P(N)      N\n\n", "");
  635.     }
  636.  
  637.     if (hdp == NULL || hdp->hitlists == NULL) {
  638.         if (hdp->hitlists == NULL) {
  639.             fprintf(fp, "\n      *** NONE ***\n");
  640.         }
  641.         return FALSE;
  642.     }
  643.  
  644. /*    dim = hdp->dim;*/ /* should be 2, or 3 for BLAST3 */
  645.  
  646.     max_id_length = get_max_ID_length(hdp->hitlists);
  647.     if (StringCmp("blastn", program) == 0 || 
  648.     StringCmp("blastp", program) == 0)
  649.         max_def_length = DEF_AND_ID_CUTOFF - max_id_length;
  650.     else
  651.         max_def_length = DEF_AND_ID_CUTOFF - max_id_length - 3;
  652.     if (max_def_length > DEF_CUTOFF)
  653.     max_def_length = DEF_CUTOFF;
  654.  
  655.     for (hlp = hdp->hitlists; hlp != NULL; hlp = hlp->next)
  656.     {
  657.         if (hlp->hsps != NULL)
  658.         {
  659.             hsp = hlp->hsps; /* first HSP */
  660.             segp = hsp->segs;
  661.             if (segp != NULL && segp->next != NULL &&
  662.                 (hit = segp->next->loc) != NULL)
  663.             {
  664.                 if (hlp->seqs != NULL)
  665.                 {
  666.                     defline = NULL;
  667.                     sdp = hlp->seqs->desc;
  668.                     if (sdp->defline != NULL) {
  669.                         defline = sdp->defline;
  670.                     }
  671.                 }
  672.  
  673.                 highScore = -1;
  674.                 TemphighScore = highScore;
  675.                 poissonScore = 99999;
  676.                 TemppoissonScore = poissonScore;
  677.                 nScore = 0;
  678.                 for (; hsp != NULL; hsp = hsp->next)
  679.                 {   /* If value associated with BLAST0_Score_sid_sum_n is 1,
  680.             then this is not set in the ASN.1 */
  681.             TempnScore = 1; 
  682.                     segp = hsp->segs;
  683.                     for (score = (ScorePtr) hsp->scores; score != NULL; score = score->next)
  684.                     {
  685.             if (score->id->str == NULL)
  686.                 scoreDescr = "";
  687.             else
  688.                 scoreDescr = score->id->str;
  689.             if (StrCmp(scoreDescr, "score") == 0)
  690.             {
  691.                             TemphighScore = score->value.intvalue;
  692.             } else if (StrCmp(scoreDescr, "p_value") == 0 ||
  693.                        StrCmp(scoreDescr, "poisson_p") == 0 ||
  694.                        StrCmp(scoreDescr, "sum_p") == 0)
  695.             {
  696.                             TemppoissonScore = score->value.realvalue;
  697.             } else if (StrCmp(scoreDescr, "poisson_n") == 0 ||
  698.                        StrCmp(scoreDescr, "sum_n") == 0)
  699.             {
  700.                             TempnScore = score->value.intvalue;
  701.             } else { /* default */
  702.             }
  703.                     }
  704.  
  705.                     if (TemphighScore > highScore)
  706.                     {
  707.                         highScore = TemphighScore;
  708.             if (segp != NULL)
  709.             {
  710.                if (StringCmp(program, "blastx") == 0)
  711.                 get_frame(segp->loc, frame_array, query_length);
  712.                else if (StringCmp(program, "tblastn") == 0)
  713.                 get_frame(segp->next->loc, frame_array, hlp->seqs->length);
  714.                else if (StringCmp(program, "tblastx") == 0)
  715.                 get_frame(segp->next->loc, frame_array, hlp->seqs->length);
  716.             }
  717.             }
  718.             if (TemppoissonScore < poissonScore)
  719.             {
  720.                         poissonScore = TemppoissonScore;
  721.                         nScore = TempnScore;
  722.                     }
  723.                 }
  724.  
  725. /* Add secondary id's (and their deflines) onto the first defline if it isn't 
  726. too long already.  */
  727.         str1 = defline;
  728.         string_length = 0;
  729.                 if (str1 != NULL)
  730.                 {
  731.                string_length = StringLen(str1);
  732.                        if (string_length <= max_def_length)
  733.             {
  734.                 if (sdp && sdp->next)
  735.                 {    
  736.                     ptr = &buffer[0];
  737.                     StrNCpy(ptr, str1, string_length);
  738.                     ptr += string_length;
  739.                     total_length = string_length;
  740.                     for (sdp1=sdp->next; sdp1; sdp1=sdp1->next)
  741.                     {
  742.                     length = format_an_id(&id[0], sdp1->id, 0);
  743.                     old_length = total_length;
  744.                     total_length += length;
  745.                     total_length += 2;
  746.                     StringCpy(ptr, " >");
  747.                     ptr += 2;
  748.                     if (total_length > max_def_length)
  749.                     {
  750. /* Add on part of last (too long) id to agree with Traditional Output.
  751. There must be some slack (extra room) in buffer for this. */
  752.                         length_left = max_def_length-old_length;
  753.                         if (length_left > 0)
  754.                         {
  755.                             StrNCpy(ptr, id, length_left);
  756.                         ptr += length_left;
  757.                         *ptr = '\0';
  758.                         }
  759.                         else
  760.                         {
  761.                             StrNCpy(ptr, "...", 3);
  762.                         ptr += 3;
  763.                         *ptr = '\0';
  764.                         }
  765.                         break;
  766.                     }
  767.                     StrNCpy(ptr, id, length);
  768.                     ptr += length;
  769.                     *ptr = ' '; ptr++;
  770.                     total_length++;
  771.                     length = StringLen(sdp1->defline);
  772.                     total_length += length; 
  773.                     if (total_length > max_def_length)
  774.                     {
  775. /* Add on part of last (too long) id to agree with Traditional Output.
  776. There must be some slack (extra room) in buffer for this. */
  777.                         length = total_length-length;
  778.                         length = max_def_length-length;
  779.                         if (length > 0)
  780.                             StrNCpy(ptr, sdp1->defline, length);
  781.                         else
  782.                             StrNCpy(ptr, sdp1->defline, 1);
  783.                         break;
  784.                     }
  785.                     StrNCpy(ptr, sdp1->defline, length);
  786.                     ptr += length;
  787.                         *ptr = '\0';
  788.                     }
  789.                     str1 = &buffer[0];
  790.                 }
  791.             }
  792.         }
  793.  
  794.                 if (str1 != NULL)
  795.                string_length = StringLen(str1);
  796.         else
  797.                string_length = 0;
  798.                 if (string_length > max_def_length)
  799.                 {
  800.                           StrNCpy (str2, str1, max_def_length);
  801.                         str2[max_def_length-3] = '.';
  802.                         str2[max_def_length-2] = '.';
  803.                         str2[max_def_length-1] = '.';
  804.                         str2[max_def_length] = '\0';
  805.         }
  806.         else
  807.         {
  808.                     if (str1 != NULL)
  809.                               StrNCpy (str2, str1, string_length);
  810.             for (index=string_length; index<max_def_length; index++)
  811.             {
  812.                 str2[index] = ' ';
  813.             }
  814.             str2[index] = '\0';
  815.         }
  816.         str1 = str2;
  817.  
  818.         
  819.  
  820.         format_an_id(&id[0], sdp->id, max_id_length);
  821.  
  822. /* Cutoff the id's after max_id_length characters to agree with Traditional 
  823. Output */
  824.  
  825.         id[max_id_length] = '\0';
  826.         if (StringCmp("blastn", program) == 0 || 
  827.             StringCmp("blastp", program) == 0)
  828.                     fprintf(fp, "%s %s %5d", id, str1, (int) highScore);
  829.         else
  830.                     fprintf(fp, "%s %s %s %5d", id, str1, frame_array, (int) highScore);
  831.  
  832. /* The following code was lifted from "print_header" in blastapp/lib/prt_hdr.c to match the traditional output. */
  833.             if (poissonScore <= 0.99) 
  834.         {
  835.                        if (poissonScore != 0.)
  836.                         fprintf(fp, "  %#-8.2lg", (double) poissonScore);
  837.                        else
  838.                         fprintf(fp, "  %#-8.2lg", 0.);
  839.             }
  840.             else if (poissonScore <= 0.999)
  841.                        fprintf(fp, "  %#-8.3lg", (double) poissonScore);
  842.             else if (poissonScore <= 0.9999)
  843.                     fprintf(fp, "  %#-8.4lg", (double) poissonScore);
  844.             else if (poissonScore <= 0.99999)
  845.                     fprintf(fp, "  %#-8.5lg", (double) poissonScore);
  846.             else if (poissonScore <  0.9999995)
  847.                     fprintf(fp, "  %#-8.6lg", (double) poissonScore);
  848.             else
  849.                     fprintf(fp, "  %#-8.7lg", (double) 1.0);
  850.         fprintf(fp, "%3d\n", (int) nScore);
  851.             }
  852.         }
  853.     }
  854.     return TRUE;
  855. }
  856.  
  857. /* The traditional BLAST output uses fasta style id's and not gi's.  Look 
  858. for a textid (i.e., fasta) and use this, otherwise use the first id. 
  859. This function operates in two modes: if id is NULL, then the length of 
  860. the id is returned; if id is not NULL, then the formatted id is returned in 
  861. "id", which should already contain space for this operation.  The ValNodePtr 
  862. vnp is actually BLAST0SeqDesc.id 
  863. The parameter max_id_length should be greater than zero if id is not NULL
  864. and one wishes to pad the end of the id with white spaces.
  865. */
  866.  
  867. static Int2
  868. format_an_id (CharPtr id, ValNodePtr vnp, Int2 max_id_length)
  869.  
  870. {
  871.     Boolean found_id;
  872.     Char temp[100];
  873.     Int2 index, length;
  874.     ValNodePtr vnp1;
  875.  
  876.     if (id != NULL)
  877.         id[0] = '\0';
  878.  
  879.     found_id = FALSE;
  880.  
  881.     for (vnp1=vnp; vnp1 != NULL; vnp1=vnp1->next)
  882.     {
  883.         if (vnp1->choice == BLAST0SeqId_textid)
  884.         {
  885.             if (id != NULL)
  886.                   sprintf(temp, "%s", vnp1->data.ptrvalue);
  887.             length = StringLen(vnp1->data.ptrvalue);
  888.             found_id = TRUE;
  889.             break;
  890.         }
  891.     }
  892.  
  893.     if (found_id == FALSE)
  894.         switch (vnp->choice) 
  895.         {
  896.         default:
  897.                 if (id != NULL)
  898.                 sprintf(temp, "Unknown");
  899.             length = 7;
  900.             break;
  901.         case BLAST0SeqId_giid:
  902.                 if (id != NULL)
  903.                 sprintf(temp, "gi|%ld", (long) vnp->data.intvalue);
  904.             length = 9;
  905.             break;
  906.         case BLAST0SeqId_textid:
  907.                 if (id != NULL)
  908.                       sprintf(temp, "%s", vnp->data.ptrvalue);
  909.                 length = StringLen(vnp1->data.ptrvalue);
  910.             break;
  911.         }
  912.  
  913.     if (max_id_length > 0)
  914.     {
  915.         if (max_id_length < length)
  916.             length = max_id_length;
  917.     }
  918.     else
  919.     {
  920.         if (BLASTCLI_ID_MAX < length)
  921.             length = BLASTCLI_ID_MAX;
  922.     }
  923.  
  924.     if (id)
  925.     {
  926.         for (index=0; index<length; index++)
  927.             id[index] = temp[index];
  928. /* Pad out the rest of id so white spaces will keep everything neat in the
  929. one-line def's */
  930.         while (index < max_id_length)
  931.         {
  932.             id[index] = ' ';
  933.             index++;
  934.         }
  935.         id[index] = '\0';    
  936.     }
  937.     return length;
  938. }
  939.  
  940. /***************************************************************************
  941. *    Get the length of an id, truncate if necessary.
  942. ***************************************************************************/
  943. static Int2 
  944. get_max_ID_length(BLAST0HitListPtr hlp)
  945.  
  946. {
  947.     BLAST0SeqDescPtr sdp;
  948.     Int2 length, id_length=0;
  949.  
  950.     for (; hlp != NULL; hlp=hlp->next)
  951.     {
  952.         if (hlp->seqs != NULL)
  953.         {
  954.             sdp = hlp->seqs->desc;
  955.             if (sdp)
  956.             {
  957.                 length = format_an_id(NULL, sdp->id, 0);
  958.                 if (length > BLASTCLI_ID_MAX)
  959.                     id_length = BLASTCLI_ID_MAX;
  960.                 if (length > id_length)
  961.                     id_length = length;
  962.             }
  963.         }
  964.     }
  965.     return id_length;
  966. }
  967.  
  968. static Int4 get_pos(BLAST0SeqIntervalPtr loc, Int4 interval)
  969. {
  970.  
  971.     Int4 position=0;
  972.  
  973.     position = loc->from + interval + 1; /* plus, default */
  974.  
  975.     switch (loc->strand)
  976.     {
  977.         case BLAST0_Seq_interval_strand_minus:
  978.         case BLAST0_Seq_interval_strand_minus_rf:
  979.             position = loc->to - interval + 1; /* minus */
  980.         break;
  981.     default:
  982.         break;
  983.     }
  984.  
  985.     return position;
  986. }
  987.  
  988. /* 
  989. "get_frame" gets the frame and puts it into an array, which the user
  990. should provide.  The array should have room for three characters.
  991. */
  992.  
  993. static void get_frame(BLAST0SeqIntervalPtr loc, CharPtr array, Int4 total_length)
  994. {
  995.  
  996.     Int2 frame;
  997.     Int4 from;
  998.  
  999.     switch (loc->strand)
  1000.     {
  1001.         case BLAST0_Seq_interval_strand_minus:
  1002.         case BLAST0_Seq_interval_strand_minus_rf:
  1003.         from = total_length-(loc->to + 1);
  1004.             frame = (Int2) (from%3) + 1;
  1005.         array[0] = '-';
  1006.         if (frame == 1)
  1007.             array[1] = '1';
  1008.         else if (frame == 2)
  1009.             array[1] = '2';
  1010.         else if (frame == 3)
  1011.             array[1] = '3';
  1012.         break;
  1013.     default:
  1014.             from = loc->from;
  1015.             frame = (Int2) (from%3) + 1;
  1016.         array[0] = '+';
  1017.         if (frame == 1)
  1018.             array[1] = '1';
  1019.         else if (frame == 2)
  1020.             array[1] = '2';
  1021.         else if (frame == 3)
  1022.             array[1] = '3';
  1023.         break;
  1024.     }
  1025.     array[2] = '\0';
  1026. }
  1027.  
  1028. /************************************************************************
  1029. *    Set up a SeqPort to allow the user to retrive residues/basepairs
  1030. *    one at a time.  Note that the BioseqPtr is really just there
  1031. *    to allow the opening of the SeqPortPtr, and is just a skeletal
  1032. *    Bioseq.
  1033. *********************************************************************/
  1034. static SeqPortPtr portnew(BioseqPtr bsp, ValNodePtr str, Int4 len)
  1035. {
  1036.     Uint1 code;
  1037.  
  1038.     bsp->length = len;
  1039.     switch (str->choice) {
  1040.         case BLAST0SeqData_ncbistdaa: 
  1041.             bsp->seq_data_type = Seq_code_ncbistdaa;
  1042.             bsp->mol = Seq_mol_aa;
  1043.             code = Seq_code_ncbieaa;
  1044.             break;
  1045.         case BLAST0SeqData_ncbi2na:
  1046.             fprintf(stderr, "Not implemented\n");
  1047.             return NULL;
  1048.         case BLAST0SeqData_ncbi4na: 
  1049.             bsp->seq_data_type = Seq_code_ncbi4na;
  1050.             bsp->mol = Seq_mol_na;
  1051.             code = Seq_code_iupacna;
  1052.             break;
  1053.         default:
  1054.             return NULL;
  1055.     }
  1056.     bsp->seq_data = str->data.ptrvalue;
  1057.     bsp->repr = Seq_repr_raw;
  1058.     return SeqPortNew(bsp, 0, -1, 0, code);
  1059. }
  1060.  
  1061. static void
  1062. TraditionalBlastWarning(BLAST0ResponsePtr blresp)
  1063.  
  1064. {
  1065.     BLAST0ResponsePtr var;
  1066.     BLAST0WarningPtr bwp;
  1067.     CharPtr warning;
  1068.  
  1069.     for (var=blresp; var; var=var->next)
  1070.     {
  1071.         if (var->choice == BLAST0Response_warning)
  1072.         {
  1073.             bwp = var->data.ptrvalue;
  1074.             warning = bwp->reason;
  1075.             wrap(stderr, "WARNING: ", warning, StringLen(warning), 79, 8);
  1076.         }
  1077.     }
  1078. }
  1079.  
  1080. #define BLSTOUT_PRINT_LEN 60
  1081.  
  1082. /************************************************************************
  1083. *    
  1084. *    Print out the alignment portion of the traditional BLAST output,
  1085. *    that is:
  1086. *
  1087. *>dbj|D32197|RURMTCB27 Mitochondrion Rupicapra rupicapra gene for cytochrome b.
  1088. *            >emb|D32197|MIRRCB27 Mitochondrion Rupicapra rupicapra gene for
  1089. *            cytochrome b.
  1090. *            Length = 646
  1091. *
  1092. *  Plus Strand HSPs:
  1093. *
  1094. * Score = 3230 (892.5 bits), Expect = 7.0e-264, P = 7.0e-264
  1095. * Identities = 646/646 (100%), Positives = 646/646 (100%), Strand = Plus / Plus
  1096. *
  1097. *Query:     1 AATACACTACACATCCGATACAGCAACAGCATTCTCCTCTGTAACCCACATTTGCCGAGA 60
  1098. *             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
  1099. *Sbjct:     1 AATACACTACACATCCGATACAGCAACAGCATTCTCCTCTGTAACCCACATTTGCCGAGA 60
  1100. *
  1101. *Query:    61 TGTAAACTACGGCTGAATCATCCGATACATACATGCAAATGGAGCATCAATATTTTTCAT 120
  1102. *             ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
  1103. *Sbjct:    61 TGTAAACTACGGCTGAATCATCCGATACATACATGCAAATGGAGCATCAATATTTTTCAT 120
  1104. *
  1105. * and so on.
  1106. ***************************************************************************/
  1107. static Boolean
  1108. TraditionalBottomBlastOutput(BLAST0ResultPtr hdp, Int4 query_length, CharPtr program, FILE *fp)
  1109. {
  1110.     SeqPortPtr spps, sppq, spp;
  1111.     Int2 num_open, num, sum_n;
  1112.     static Bioseq bios, bioq;
  1113.     BLAST0SegmentPtr seg;
  1114.     BLAST0SeqIntervalPtr loc, loc1, loc2, loc_next, loc_s;
  1115.     BLAST0HitListPtr hlp;
  1116.     Int4 nbr_ident, nbr_pos, pos_temp;
  1117.     BLAST0SeqDescPtr sdp;
  1118.     Char ch, chh, qchh, chh_std, qchh_std;
  1119.     Char frame_array[3];
  1120.     CharPtr defline, cp, head, s;
  1121.     CharPtr strand1, strand2;
  1122.     ScorePtr score;
  1123.     CharPtr scoreDescr;
  1124.     BLAST0HSPPtr hsp;
  1125.     int i, sc=0, ldef;    /* Left as "int" for use in Warren's wrap */
  1126.     FloatHi pv=0, ev=0;
  1127.     Int4 j, pos, range, spos, from, to; 
  1128.     Int4 first, width;
  1129.     Int4 hsp_length, qposition, sposition;
  1130.     Uint2 prev_strand;
  1131.     ValNodePtr vnp;
  1132.     
  1133.     
  1134. /* Start of main loop that runs throughout the function. */
  1135.     for (hlp = hdp->hitlists; hlp != NULL; hlp = hlp->next)     {
  1136.         if ((hsp = hlp->hsps) == NULL) {
  1137.             continue;
  1138.         }
  1139.         if (hlp->seqs == NULL) {
  1140.             continue;
  1141.         }
  1142.         if ((sdp = hlp->seqs->desc) == NULL) {
  1143.             continue;
  1144.         }
  1145.         ldef = 0;
  1146.         for (sdp = hlp->seqs->desc; sdp != NULL; sdp = sdp->next) {
  1147.             ldef += strlen(sdp->defline) + 1;
  1148.             ldef += format_an_id(NULL, sdp->id, 0);
  1149.             ldef += 2;
  1150.         }
  1151.         s = defline = MemNew(ldef + 1);
  1152.         for (sdp = hlp->seqs->desc; sdp != NULL; sdp = sdp->next) 
  1153.         {
  1154.             if (s != defline) {
  1155.                 sprintf(s++, " ");
  1156.             }
  1157.             sprintf(s++, "%c", '>');
  1158. /* Don't pad the id here, as it's secondary and goes into the def line */
  1159.             format_an_id(s, sdp->id, 0);
  1160.             s += StringLen(s);
  1161.             sprintf(s, " %s", sdp->defline);
  1162.             s += StringLen(s);
  1163.         }
  1164.         
  1165.         if (defline != NULL) {
  1166.             for (cp = defline; ((ch = *cp) != NULLB && !IS_WHITESP(ch)); cp++);
  1167.             for (; IS_WHITESP(*cp); cp++) ;
  1168.             i = cp - defline;
  1169.             i = MIN(i, 12);
  1170.             fprintf(fp, "\n\n");
  1171.             wrap(fp, "", defline, ldef, 79, i);
  1172.             defline = MemFree(defline);
  1173.         }
  1174.         fprintf(fp, "%*sLength = %s\n", i, "", 
  1175.             Ltostr(hlp->seqs->length, 1));
  1176.         prev_strand = UINT2_MAX;    /* Value that will not be used*/
  1177.         for (; hsp != NULL; hsp = hsp->next) {
  1178.             sum_n = 1;    /* Default (one) not in structure */
  1179.             for (score = (ScorePtr) hsp->scores; score != NULL; score = score->next) {
  1180.                 if (score->id->str == NULL)
  1181.                     scoreDescr = "";
  1182.                 else
  1183.                     scoreDescr = score->id->str;
  1184.                 if (StrCmp(scoreDescr, "score") == 0)
  1185.                 {
  1186.                                 sc = score->value.intvalue;
  1187.                 } else if (StrCmp(scoreDescr, "p_value") == 0 ||
  1188.                            StrCmp(scoreDescr, "poisson_p") == 0 ||
  1189.                            StrCmp(scoreDescr, "sum_p") == 0)
  1190.                 {
  1191.                                 pv = score->value.realvalue;
  1192.                 } else if (StrCmp(scoreDescr, "e_value") == 0 ||
  1193.                            StrCmp(scoreDescr, "poisson_e") == 0 ||
  1194.                            StrCmp(scoreDescr, "sum_e") == 0)
  1195.                 {
  1196.                                 ev = score->value.realvalue;
  1197.                 } else if (StrCmp(scoreDescr, "poisson_n") == 0 ||
  1198.                            StrCmp(scoreDescr, "sum_n") == 0)
  1199.                 {
  1200.                                 sum_n = score->value.intvalue;
  1201.                 } else { /* default */
  1202.                 }
  1203.             }
  1204.             hsp_length = hsp->len;
  1205. /* loc1 is the query; loc2 is the subject. */
  1206.             loc1 = hsp->segs->loc;
  1207.             loc2 = hsp->segs->next->loc;
  1208.             if (StringCmp(program, "blastn") == 0 ||
  1209.                 StringCmp(program, "blastp") == 0)
  1210.             {
  1211.                 loc = loc1;
  1212.                 loc_s = loc2;
  1213.             }
  1214.             else if (StringCmp(program, "blastx") == 0) 
  1215.             {
  1216.                 loc = loc2;
  1217.                 loc_s = loc1;
  1218.             }
  1219.             else if (StringCmp(program, "tblastn") == 0) 
  1220.             {
  1221.                 loc = loc1;
  1222.                 loc_s = loc2;
  1223.             }
  1224.             else if (StringCmp(program, "tblastx") == 0) 
  1225.             {
  1226.                 loc = loc1;
  1227.                 loc_s = loc2;
  1228.             }
  1229.  
  1230.             if (loc_s->strand != 0 && loc_s->strand != prev_strand) 
  1231.             {
  1232.                 prev_strand = loc_s->strand;
  1233.                 fprintf(fp, "\n  %s Strand HSPs:\n", 
  1234.                     gish_str_strand(loc_s->strand));
  1235.             }
  1236.             fprintf(fp, "\n Score = %ld (%.1lf bits), Expect = %#0.2lg, ",
  1237.                 (long)sc, (double) sc * la2, (double) ev);
  1238.  
  1239.             if (sum_n < 2) {
  1240.                 fprintf(fp, "P = %#0.2lg\n", (double) pv);
  1241.             } else {
  1242.                 fprintf(fp, "Sum P(%u) = %#0.2lg\n", sum_n, (double) pv);
  1243.             }
  1244.  
  1245.             nbr_ident = nbr_pos = 0;
  1246.             if ((sppq = portnew(&bioq, hsp->segs->str, hsp_length)) == NULL) {
  1247.                 return FALSE;
  1248.             }
  1249.             if ((spps = portnew(&bios, hsp->segs->next->str, hsp_length)) == NULL) {
  1250.                 return FALSE;
  1251.             }
  1252.             
  1253.             range = loc->to - loc->from;
  1254.             if (StringCmp(program, "tblastx") == 0) 
  1255.                 range /= 3;
  1256.             for (pos=0; pos <= range; pos++) 
  1257.             {
  1258.                 chh = SeqPortGetResidue(spps);
  1259.                 qchh = SeqPortGetResidue(sppq);
  1260.                 if (chh == qchh) 
  1261.                 {
  1262.                     nbr_ident++;
  1263.                     if (StringCmp(program, "blastn") == 0) 
  1264.                         nbr_pos++;
  1265.                 }
  1266.                 if (StringCmp(program, "blastn") != 0) 
  1267.                 {
  1268.                      chh_std = ralpha[ chh];
  1269.                     qchh_std = ralpha[qchh];
  1270.                     if (matr[chh_std][qchh_std] > 0) 
  1271.                         nbr_pos++;
  1272.                 }
  1273.             }
  1274.  
  1275.             fprintf(fp,
  1276.             " Identities = %ld/%ld (%ld%%), Positives = %ld/%ld (%ld%%)",
  1277.             (long) nbr_ident, (long) pos, (long)((100*nbr_ident)/pos),
  1278.             (long) nbr_pos, (long) pos, (long)((100*nbr_pos)/pos) );
  1279.             if (StringCmp( program, "blastn") == 0)
  1280.             {
  1281.                 strand1 = gish_str_strand(loc1->strand);
  1282.                 strand2 = gish_str_strand(loc2->strand);
  1283.                 if (StringCmp(strand1, "Undefined") != 0 &&
  1284.                     StringCmp(strand2, "Undefined") != 0)
  1285.                     fprintf(fp, ", Strand = %s / %s", strand1, strand2);
  1286.             }
  1287.             else if (StringCmp( program, "blastx") == 0)
  1288.             {
  1289.                 loc = hsp->segs->loc;
  1290.                 get_frame(loc, frame_array, query_length);
  1291.                 fprintf(fp, ", Frame = %s", frame_array);
  1292.             }
  1293.             else if (StringCmp( program, "tblastn") == 0)
  1294.             {
  1295.                 loc = hsp->segs->next->loc;
  1296.                 get_frame(loc, frame_array, hlp->seqs->length);
  1297.                 fprintf(fp, ", Frame = %s", frame_array);
  1298.             }
  1299.             else if (StringCmp( program, "tblastx") == 0)
  1300.             {
  1301.                 loc = hsp->segs->loc;
  1302.                 get_frame(loc, frame_array, query_length);
  1303.                 fprintf(fp, ", Frame = %s", frame_array);
  1304.                 loc = hsp->segs->next->loc;
  1305.                 get_frame(loc, frame_array, hlp->seqs->length);
  1306.                 fprintf(fp, " / %s", frame_array);
  1307.             }
  1308.             putc('\n', fp);
  1309.             num_open = 1;
  1310. /* Print out the actual alignment. */ 
  1311.             for (pos=0; pos < hsp_length; pos += BLSTOUT_PRINT_LEN) {
  1312.                 head = "Query";
  1313.                 first = TRUE;
  1314.                 num = 0;
  1315.                 for (seg = hsp->segs; seg != NULL; seg = seg->next, num++) {
  1316.                     loc = seg->loc;
  1317.                     from = loc->from;
  1318.                     to = loc->to;
  1319.                     if (first)
  1320.                     {
  1321.                          qposition = get_pos(loc, pos);
  1322.                          width = MAX(5, Lwidth(qposition, 1));
  1323.                          loc_next = seg->next->loc;
  1324.                          sposition = get_pos(loc_next, pos);
  1325.                          width = MAX(width, Lwidth(sposition, 1));
  1326.                     }
  1327.                     if (!first) {
  1328.                         SeqPortSeek(sppq, pos, 0);
  1329.                         if (num_open != num) {
  1330.                             spps = portnew(&bios, seg->str, hsp_length);
  1331.                             if (spps == NULL) {
  1332.                                 return FALSE;
  1333.                             }
  1334.                             num_open = num;
  1335.                         }
  1336.                         SeqPortSeek(spps, pos, 0);
  1337.                         for (j=0; j<(width+8); j++)
  1338.                             putc(' ', fp);
  1339. /*
  1340.                         fprintf(fp, "%s", " ");
  1341. */
  1342.                         for (j=0; j+pos < hsp_length && j < BLSTOUT_PRINT_LEN; j++) {
  1343.                              chh = SeqPortGetResidue(spps);
  1344.                             qchh = SeqPortGetResidue(sppq);
  1345.                             if (prots) {
  1346.                                 if (chh == qchh) {
  1347.                                     putc(chh, fp);
  1348.                                 } else {
  1349.                                     chh_std = ralpha[ chh];
  1350.                                     qchh_std = ralpha[qchh];
  1351.                                     putc((matr[chh_std][qchh_std]>0) ? '+' : ' ', fp);
  1352.                                 }
  1353.                             } else {
  1354.                                 putc((chh == qchh) ? '|' : ' ', fp);
  1355.                             }
  1356.                         }
  1357.                         spp = spps;
  1358.                     } else {
  1359.                         spp = sppq;
  1360.                     }
  1361.                     fprintf(fp, "\n");
  1362.                     SeqPortSeek(spp, pos, 0);
  1363.                     if (pos+BLSTOUT_PRINT_LEN < hsp_length)
  1364.                         spos = BLSTOUT_PRINT_LEN - 1;
  1365.                     else
  1366.                         spos = hsp_length - pos - 1;
  1367.                     pos_temp = pos;
  1368.                         if (StringCmp( program, "tblastx") == 0)
  1369.                     {
  1370.                         pos_temp *= 3;
  1371.                     }
  1372.                     else if (StringCmp(head, "Query") == 0)
  1373.                     {
  1374.                          if (StringCmp(program, "blastx") == 0)
  1375.                          {
  1376.                             pos_temp *= 3;
  1377.                          }
  1378.                     }
  1379.                     else if (StringCmp(head, "Sbjct") == 0)
  1380.                     {
  1381.                          if (StringCmp(program, "tblastn") == 0)
  1382.                          {
  1383.                             pos_temp *= 3;
  1384.                          }
  1385.                     }
  1386.                     qposition = get_pos(loc, pos_temp);
  1387.                     fprintf(fp, "%s: %*ld ", head, (int) width, (long) qposition);
  1388.                     for (j=0; j <= spos; j++) {
  1389.                         chh = SeqPortGetResidue(spp);
  1390.                         fputc(chh ,fp);
  1391.                     }
  1392.                     pos_temp = j-1+pos;
  1393.                         if (StringCmp( program, "tblastx") == 0)
  1394.                     {
  1395.                         pos_temp *= 3;
  1396.                         pos_temp += 2;
  1397.                     }
  1398.                     else if (StringCmp(head, "Query") == 0)
  1399.                     {
  1400.                          if (StringCmp(program, "blastx") == 0 ||
  1401.                                   StringCmp( program, "tblastx") == 0)
  1402.                         {
  1403.                             pos_temp *= 3;
  1404.                             pos_temp += 2;
  1405.                         }
  1406.                     }
  1407.                     else if (StringCmp(head, "Sbjct") == 0)
  1408.                     {
  1409.                          if (StringCmp(program, "tblastn") == 0)
  1410.                          {
  1411.                             pos_temp *= 3;
  1412.                             pos_temp += 2;
  1413.                          }
  1414.                     }
  1415.                     qposition = get_pos(loc, pos_temp);
  1416.                     fprintf(fp, " %ld\n", (long) qposition);
  1417.                     head = "Sbjct";
  1418.                     first = FALSE;
  1419.                 }
  1420.             }
  1421.             SeqPortFree(sppq);
  1422.             SeqPortFree(spps);
  1423.         }
  1424.     }
  1425.     fprintf(fp, "\n");
  1426.     return TRUE;
  1427. }
  1428.  
  1429. /********************************************************************
  1430. *    
  1431. *    Formats the parameters (e.g., command line options) sent 
  1432. *    back by the server and the statistics found at the bottom
  1433. *    of the output.
  1434. *********************************************************************/
  1435. static Boolean
  1436. TraditionalTailBlastOutput(BLAST0ResponsePtr blresp, FILE *fp)
  1437. {
  1438.     ValNodePtr res;
  1439.     ValNodePtr parms, stats;
  1440.     CharPtr s;
  1441.     Boolean new;
  1442.  
  1443.     if ((res = find(blresp, BLAST0Response_parms)) != NULL) {
  1444.         if ((parms = res) != NULL) {
  1445.             fprintf(fp, "\nParameters:\n");
  1446.             new = TRUE;
  1447.             for (; parms != NULL; parms = parms->next) {
  1448.                 s = parms->data.ptrvalue;
  1449.                 if (*s == '\0') {
  1450.                     fprintf(fp, "\n", s);
  1451.                     new = TRUE;
  1452.                 } else {
  1453.                     if (new) {
  1454.                         fprintf(fp, "  ");
  1455.                     }
  1456.                     fprintf(fp, "%s", s);
  1457.                     new = FALSE;
  1458.                 }
  1459.             }
  1460.         } else {
  1461.             fprintf(fp, "No parameters received from server\n");
  1462.         }
  1463.     }
  1464.     if ((res = find(blresp, BLAST0Response_stats)) != NULL) {
  1465.         if ((stats = res) != NULL) {
  1466.             fprintf(fp, "\nStatistics:\n");
  1467.             new = TRUE;
  1468.             for (; stats != NULL; stats = stats->next) {
  1469.                 s = stats->data.ptrvalue;
  1470.                 if (*s == '\0') {
  1471.                     fprintf(fp, "\n", s);
  1472.                     new = TRUE;
  1473.                 } else {
  1474.                     if (new) {
  1475.                         fprintf(fp, "  ");
  1476.                     }
  1477.                     fprintf(fp, "%s", s);
  1478.                     new = FALSE;
  1479.                 }
  1480.             }
  1481.         } else {
  1482.             fprintf(fp, "No statistics received from server\n");
  1483.         }
  1484.     }
  1485.     return TRUE;
  1486. }
  1487.  
  1488. Boolean LIBCALL TraditionalBlastOutput(BLAST0ResultPtr hdp, BLAST0ResponsePtr blresp, CharPtr program, FILE *fp)
  1489. {
  1490.     Int4 query_length=0;
  1491.  
  1492.     if (!TraditionalHeadBlastOutput(blresp, &query_length, program, fp)) {
  1493.         return FALSE;
  1494.     }
  1495.    
  1496.    if (!TraditionalHistBlastOutput(hdp, fp)) {
  1497.       return FALSE;
  1498.    }
  1499.    
  1500.   
  1501.     if (TraditionalTopBlastOutput(hdp, query_length, program, fp)) 
  1502.     {
  1503.         TraditionalBlastWarning(blresp);
  1504.         if (!TraditionalBottomBlastOutput(hdp, query_length, program, fp)) {
  1505.             return FALSE;
  1506.         }
  1507.     }
  1508.     if (!TraditionalTailBlastOutput(blresp, fp)) {
  1509.         return FALSE;
  1510.     }
  1511.     return TRUE;
  1512. }
  1513.